Role of anisotropy for protein-protein encounter 
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Protein-protein interactions comprise both transport and reaction steps. During the transport 
step, anisotropy of proteins and their complexes is important both for hydrodynamic diffusion 
and accessibility of the binding site. Using a Brownian dynamics approach and extensive computer 
simulations, we quantify the effect of anisotropy on the encounter rate of ellipsoidal particles covered 
with spherical encounter patches. We show that the encounter rate k depends on the aspect ratios 
£ mainly through steric effects, while anisotropic diffusion has only a little effect. Calculating 
analytically the crossover times from anisotropic to isotropic diffusion in three dimensions, we find 
that they are much smaller than typical protein encounter times, in agreement with our numerical 
results. 
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Protein-protein interactions are at the heart of most 
molecular processes in biological systems and are inten- 
sively investigated both by experiment and theory [lj. 
Conceptually, protein-protein binding consists of trans- 
port and reaction steps. Brownian dynamics simulations 
have been introduced to study transport towards asso- 
ciation of spherical model particles @, Isotropy of 
the diffusion process can be used to develop computer 
time efficient methods for propagating reactive particles 
in time and space However, assuming spherical par- 
ticles does not address two important aspects of protein 
encounter. First, proteins and their complexes are typ- 
ically not spherical, but their shape and diffusional be- 
haviour can be highly anisotropic. Second, their bind- 
ing sites are strongly localized and thus binding is not 
isotropic neither. Both of these issues have been ad- 
dressed in Brownian dynamics simulations, e.g. by us- 
ing approximate schemes for the mobility tensor of arbi- 
trarily shaped proteins [5[ and by incorporating all-atom 
force fields However, the full role of anisotropy for 
protein encounter has not been systematically studied be- 
fore in a generic model for protein-protein encounter. Al- 
though not suited to answer specific biological questions, 
simple models without atomic details also have the ad- 
vantage that they can be easily upscaled to large system 
sizes. 

In order to study theoretically the effect of shape 
anisotropy on protein-protein encounter, we systemati- 
cally vary the aspect ratio £ = R»/R± of ellipsoids with 
a rotational symmetry, compare Fig.QJa). This is one of 
the few systems for which closed analytic expressions of 
the translational (') and rotational ( r ) friction coefficients 
d/, r ± (R\\,Rx) are known [tHH- Recently, the diffusional 
properties of such anisotropic particles have been inves- 
tigated in great detail for two dimensions [t| [l(J ■ Here 
we study three dimensions, both numerically and ana- 
lytically. In addition, we not only consider anisotropy 
in diffusion, but also anisotropy in binding by placing 
spherical encounter patches at various positions on the 
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FIG. 1: (a) Geometry of an ellipsoidal model particle 
with a reactive patch on its surface at some angle tp 
with respect to the symmetry axis eii . (b) Illustration of 
an encounter configuration. 



ellipsoids. This approach allows us to assess in quanti- 
tative detail the relative role of hydrodynamic and steric 
anisotropy for protein-protein encounter. 

To quantify the effect of anisotropy on molecular en- 
counter, we performed Brownian dynamics simulations 
with a pair of ellipsoids in a periodic boundary box of 
edge length L, which represents the concentration of our 
system. The appropriate Langevin equation is evolved in 
discrete time steps At via the Euler algorithm [Tlj : 



X(t + At) = X(t) + g(At) + 0{At 2 ) 



(1) 



All vectors are six-dimensional generalized coordinates 
including the orientational information. The noise term 
g(Af) is determined by the Einstein relation (g) = 0, 
(g(At)-g t (At)) = 2k B TJAAt, where T a is ambient tem- 
perature and M = <^ _1 is the 6x6 mobility matrix. In 
the following, all lengths are given in units of Ur, = R± 
and times are scaled by Ut — C, r [R^_, R\_) / {k B T a ). In 
the simulations we chose Ul — 2nm and thus Ut ~ 50ns 
at viscosity of water 77 = ImPas. As time step we use 
At = 2 • 1CP 5 . We do not consider direct forces between 
the model particles in this study; in particular, we ne- 
glect electrostatic and two-body hydrodynamic interac- 
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FIG. 2: (a)-(c) Encounter rates for different patch sizes and locations, (d)-(f) Statistical estimate of the fraction of 
configuration space available for encounter. The patches are located according to 01 in (a) and (d), 02 in (b) and 

(e) and 03 in (c) and (f). 



tions. Otherwise we implement hard-core repulsion. An 
analytic overlap criterion for a pair of ellipsoids is difficult 
to derive, but suitable algorithms based on the solution 
of the characteristic equation of the two ellipsoids have 
been derived for hard body fluids 

Each of the model particles carries a spherical reac- 
tive patch of radius r on its surface whose location is 
described by the angle ip, compare Fig.QJa). An overlap 
of these reactive patches is considered as an encounter, 
compare Fig.QJb). We average over all initial conditions 
by starting 10 simulations at random initial positions 
and orientations for each parameter set. The main quan- 
tity of interest is the encounter rate k, defined as the in- 
verse first passage time to encounter. As it is common in 
Brownian dynamics of proteins, the encounter rate thus 
emerges from the definition of an absorbing boundary for 
the random walk For the two model particles in the 
simulation box we consider three scenarios of patch lo- 
cations: fa = f 2 = (01); ip! = 0, tp 2 = tt/2 (02); 

In a coordinate system spanned by the principal axes, 
the friction matrix of an ellipsoid is diagonal. Taking the 

corresponding values for a sphere of radius R± as a refer- 

t It 

ence scale, the diffusion coefficients D,/, ± only depend on 
£. The relative translational mobility of an ellipsoid com- 
pared to a sphere is then given as A(£) = Dt (£) + 2.Dj_(£). 
According to the Smoluchowski equation, the encounter 
rate is expected to depend linearly on the mobility and 
the concentration of the reacting particles. Therefore we 



normalize the encounter rates obtained by our simula- 
tions by A(£) and by the volume of our simulation box 
(here L 3 = 100 3 ). The corresponding data is shown in 
Figs. [2Ja)~(c). The rates are given in 1/Ut- In each 
case, we exclude aspect ratios for which the reaction 
patches span the whole ellipsoid. Our simulation re- 
sults show that encounter rates can vary up to two or- 
ders of magnitude depending on aspect ratio and patch 
position. While for patches located at the tips (01), en- 
counter efficiency increases with aspect ratio, it decreases 
for patches located at the sides (03). For the mixed case 
(02), changes in aspect ratio have only a weak effect. 

When interpreting these results, an important system- 
atic difference between different aspect ratios £ has to 
be noted: As the geometry of the ellipsoid is changing 
with £, the exposed volume fraction fy of the reactive 
patches not covered by the steric particle is changing. 
We estimated the steric effect on the encounter rate by 
studying the fraction of non-overlapping ellipsoid config- 
urations with touching reactive patches / no . A scheme 
of the setup is shown in Fig. QJb). Particularly, the cen- 
ters of the patches are placed at the distance 2r and 
both ellipsoids are randomly rotated around the cen- 
ter of their respective patch. The results for / no from 
drawing 10 5 of such random encounter configurations are 
shown in Figs.[2]^d)-(f) for the parameters corresponding 
to Figs.[2ja)-(c). The plots show that the qualitative fea- 
tures of the encounter rates are well reproduced by the 
steric constraints. This leads to the conclusion that the 
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main reason for the changes in the encounter rate in the 
preceeding study is not the altered hydrodynamic behav- 
ior of the ellipsoids but the steric hindering of encounters 
due to the changing geometry. 

To obtain another measure for the effect of hydro- 
dynamics, we compared our simulations with explicit 
anisotropic diffusion to simulations, where we do not ac- 
count for the anisotropy in the diffusion matrix. Par- 
ticularly, we use an isotropic diffusion matrix with the 
average translational and rotational diffusion coefficients 
f* vg = + 2Z? t L )/3, which is equal to the isotropic 
limit, and -Davg = (DT, + 2D r L )/3, respectively. We per- 
form the same simulations as in Fig. [2] with this new D. 
Moreover, we compared to simulations at higher effective 
concentration (L = 20). Fig. [3] shows the relative devia- 
tion of the encounter rates A iso = fei so tropicAanisotropic— 1- 
Interestingly, in case <?3 there is no significant deviation 
from the original results. However, considering Q\ the 
(artificial) isotropic encounter rate is larger for prolates 
(£ > I). This effect can also be observed in the mixed 
case Q2, decreased by roughly a factor 2. This is reason- 
able as here only one of the two encountering ellipsoids 
has its patch at %j) = 0. These findings show that the 
anisotropic diffusion of elongated ellipsoids leads to a de- 
crease of the encounter rates. However, the deviation up 
to £ = 5 is moderate (A; so w 0.7), so we conclude that the 
impact of anisotropic diffusion on molecular encounter is 
rather weak. 

Anisotropic diffusion in isotropic environments is only 
relevant on small time and length scales. In the follow- 
ing, we derive the effective translational diffusion prop- 
erties for a finite time window At considering an arbi- 
trary body with three different translational (Dj) and 
rotational (Df) diffusion coefficients along the principal 
axis (in the body-fixed coordinate system). As the ro- 
tations due to rotational diffusion are supposed to be 
completely independent, the crossover only affects the 
effective translational diffusion which is described by a 
3x3 diffusion matrix: D' - = SijDi, where Sij is the 
Kronecker delta. A rotation determined by a vector of 
angles 6 = {9i,9 2 ,9 3 ) J ' around the three principal axes 
can be described by the rotation matrix S = e~ ej , with 
J = (J 1 , J 2 , J 3 )^. 9 ■ J denotes a formal scalar product 
and I k are matrices defined by Jy = e,jy , where e is the 
Levi-Civita symbol. We proceed considering only small 
rotations occurring at small times, so that we can expand 
the rotation in orders of 9: 

S« = h 2 -f- + e ikj k 6 -=^ + °-f- + 0(6*) , (2) 

where 2 = Y^=i ®\- This rotation is applied to D* and 
we get W = §D*S T , where we will only consider terms 
up to second order in 9 in the following. Furthermore 
we average over all possible orientations, weighted by 
the probability density p(9i,t) ~ exp(— 0?/(4D?t)) due 
to rotational diffusion. As we consider only small 9 and 



t <C 1, it will now also be sufficient to expand p(9i,t) 
up to second order in t. Hence, the Gaussian probability 
distribution can be replaced by a uniform distribution 
regarding correct integral boundaries. The non-diagonal 
entries are odd in 9i and thus vanish. The average of the 
diagonal entries is the central quantity for the calculation 
of the mean square displacement: 

toi/2 w 2 /2 w 3 /2 

(B*u)t = / cWi / d9 2 f d9 3 W tl , (3) 

WiW 2 W 3 J J J 

-Mil/2 -1C2/2 -UJ3/2 

where Wi — \j2AD\t is the width of the uniform dis- 
tribution interval. Only terms of even orders of 9i will 
contribute. That is, the integral in Eq. [3] will only lead 
to zeroth and second moments of the angular distribu- 
tion. By the average action of the rotational diffusion, Di 
transforms into an effective translational diffusion con- 
stant (B*jj)t over time t in the fixed laboratory coordi- 
nate space. Considering a vector D = (Di,D 2 ,D 3 ), the 
evolution of the effective, orientation averaged vector of 
diffusion coefficients (D)(t) = «B*u) t , (B t 2 2 )t, (B*33>t) 
can be expressed in a matrix form, not taking into ac- 
count non-linear terms in t: 

(D)(t) = D-R(t)+0(t 2 ) , (4) 
%(t) = 2t\e ijk \D e k + <%(1 + 2t(Df - D°)) , (5) 

where D e — Y^k=i The principal axes of effective 
motion are constant since the coupling terms between 
the translational degrees of freedom vanish when aver- 
aging over all possible orientations. Because rotational 
diffusion is independent of time and orientation, R keeps 
its form for all times. We can also apply M.(St) to some 
diffusion vector (D)(t). Thus, it is possible to evaluate 
the effective change of the diffusion coefficients for large 
times t in small steps St = t/N with only making errors 
of 0(N5t 2 ): 

(D)(N6t) = D ■ (R(5t)f + 0(NSt 2 ) . (6) 

In the limit of large N the error in Eq. ([6]) vanishes: 
0(Nt 2 /N 2 ) = 0{t 2 /N) -> 0. This basically means that 
we calculate (D)(i) with infinite accuracy. We can now 
evaluate the overall, orientation averaged mean square 
displacement by integrating (D) (t) with respect to t: 

(2ft + (ad 2 — bd 3 )g- + cdi5 + \ 
2ft + (ad 1 -bd2)g- + cd 2 g+ , (7) 
2ft + {ad 3 - bd{)g- + cd 3 g+ J 

with a = Dl - Df, b = D\ - D e 2l c = Va 2 + b 2 + ab, 
d t = 3(A - /), / = Ei=i A/3, and 

1 _ „-2(Zje-e)t 1 _ -2(D°+c)t 

g± = ± • (8) 

y 6(D e -c)c 6{D e + c)c V ; 
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FIG. 3: (a) Relative deviation of the encounter rates from the original data (compare Fig. [5]) assuming isotropic 
motion at all times. Reactive patches have a radius of r = 0.5. The dotted line indicates the supposed scaling and 
the dashed line marks Ai SO = 0. Data from simulations with L = 100 is shown with red, filled symbols, data with 
L = 20 is shown with green, hollow symbols, (b) Example of the time window dependency of the effective principal 
translational diffusion coefficients for the following choice of parameters (generic units): D\ — 3, Z?2 = 2, D$ = 1, 

Df = 0, Z?2 = 0.005, D\ = 0.5. The data points have been obtained by simulation (the error bars depict the 
standard deviation obtained from 10 5 individual runs), the solid lines represent the theoretic prediction and the 
dashed lines indicate the two relevant time scales for the crossover as given in Eq. ©. (c) Crossover times and 
rotational friction coefficients for ellipsoids for a wide range of aspect ratios. 



This result can be used to obtain effective translational 



diffusion coefficients via D cff (Ai) 



'■(At))/2t, 



pare Fig. [3jb). Thus, the crossover from anisotropic 
to isotropic diffusion in 3D occurs on two time scales 
Ty 2 = 1/ (2D e ± c) which increasingly deviate with in- 
creasing anisotropy. 

The corresponding crossover times for ellipsoids are 
Tf = l/(4Dff(0 + 2£>i(0) and T 2 C = 1/(6£>1(£)). Th e 
rotational friction coefficients and their implication for 
Tw 2 are shown in Fig. Etc). The asymptotic behavior 
indicated in the plot can be derived from the full solu- 
tion of the friction coefficients. For £ < both and 
D 1 ^ approach a constant value. Therefore, regarding ro- 
tational diffusion no significant differences arc expected 
for oblates, which corresponds well to the findings from 
Fig. [3ta). In contrast, rotational diffusion particularly 
around ej_, which governs Tlf, is strongly decreased for 
large £. Therefore, the relevant range of anisotropic dif- 
fusion grows and the indication of the scaling in Fig. [3ja) 
shows that this again corresponds well to the simulation 
data. Assuming this scaling to govern the impact of hy- 
drodynamic anisotropy, one might argue that the effect 
will be strongest for £ 3> 1 and small L, i.e. large con- 
centrations. However, if we require L > 2£, so that the 
particles fit into the simulation box, encounter times are 
larger than T% for all values of £. Thus our analytical 
calculation confirms that anisotropic diffusion does not 
have a strong effect on protein-protein encounter rates. 
This finding also validates computational schemes which 
assume isotropic diffusion for efficient description of the 
diffusion steps 

An interesting question that has not been addressed 
yet is whether the effect of hydrodynamic anisotropy is 



stronger regarding re-encounter. Dynamic dissociation 
and re-association can be easily realized in our model for 
example by introducing a finite average lifetime for each 
bond. Recently it has been shown with a similar Brow- 
nian dynamics approach that productive protein-protein 
encounter is preceded by many non-productive contacts 
11]. In a similar vein, dissociation after binding is likely 
to be followed by additional encounters. Because pre- 
lates have a lower overall rotational diffusion coefficient, 
one expects that they are more likely to return to an 
encounter. It might well be that for re-encounter, the 
aspect ratio plays a more important role for accessibility 
of the binding site than found here for protein diffusion 
to the first encounter. Similar considerations might be 
valid for protein dynamics in the context of protein clus- 
ters, where diffusion might be restricted by the presence 
of other components. 
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